%!TEX TS-program = xelatex
%!TEX encoding = UTF-8 Unicode
\documentclass[12pt]{article}
\usepackage[margin=2cm]{geometry}
\geometry{a4paper}
\usepackage{graphicx}
\usepackage{amssymb}
\usepackage{fontspec,xltxtra,xunicode,amsmath}
\defaultfontfeatures{Mapping=tex-text}
\setmainfont{文泉驿正黑}
\usepackage{xeCJK}
\usepackage{indentfirst}
\setlength{\parindent}{2.22em}
\usepackage{pdfpages}
\usepackage{listings}
\lstset{
        language=C++,
        keywordstyle=\bfseries\ttfamily\color[rgb]{0,0,1},
        identifierstyle=\ttfamily,
        commentstyle=\color[rgb]{0.133,0.545,0.133},
        stringstyle=\ttfamily\color[rgb]{0.627,0.126,0.941},
        showstringspaces=false,
        basicstyle=\small,
        numberstyle=\footnotesize,
        numbers=left,
        stepnumber=1,
        numbersep=10pt,
        tabsize=2,
        breaklines=true,
        prebreak = \raisebox{0ex}[0ex][0ex]{\ensuremath{\hookleftarrow}},
        breakatwhitespace=false,
        aboveskip={1.5\baselineskip},
  columns=fixed,
  upquote=true,
  extendedchars=true,
% frame=single,
% backgroundcolor=\color{lbcolor},
}
\usepackage[hidelinks=true]{hyperref}
\newcommand{\ud}{\mathrm{d}}
\renewcommand\refname{参考文献}
\renewcommand{\contentsname}{目录}
\renewcommand{\abstractname}{摘要}
\title{计算概论：：流体世界\\ Project Fluid World}
\author{徐智怡 \and 黄康靖 \and 沈钟灵}
\date{}                                         
\setcounter{secnumdepth}{5}
\setcounter{tocdepth}{5}
\begin{document}
\maketitle
\begin{abstract}
为了进一步在实践中掌握c++语言，理解面向对象编程；应用流体力学知识并初步了解一些数值方法的相关内容。
以c++语言为编程工具，数值模拟并可视化展示钝物体绕流的过程。对于一个粘性流体绕圆柱的二维层流流动问题，采用流函数-涡量法的高阶混合有限差分格式、初始条件和边界条件的处理方法。可视化的部分则采用染色线生成方法。\end{abstract}
\tableofcontents
\newpage
\section{课题目的}
我们常常惊叹于流体的神秘莫测，它和生活本身一样难以预测。作为经典力学的最后一块堡垒，科学家在上个世纪下半叶找到了一种极为有效的攻克它的途径——计算流体力学\footnote{即computational fluid dynamics,CFD}。随着计算机科学的不断发展，更强大的硬件与更巧妙的算法让这门方法逐渐成熟，能够解决许多工业上的问题。介于这门学科的复杂精深，我们想要实现的功能是单一的。即便如此，这个课题对我们的挑战，起到的提高促进作用仍是巨大的。

在本课题中，我们要实现的目标是：
\begin{enumerate}
\item 采用SIMPLE算流函数涡量高阶混合有限差分格式数值模拟二维层流流动。
\item 利用openGL，Qt提供的接口实现计算结果的可视化。
\item 由于计算时间漫长。在程序主要结构中镶嵌一个数独游戏，以供使用者消遣。
\end{enumerate}

\section{实现方案}
\subsection{流体力学基本原理}
流体作为连续介质，在宏观尺度分析时可以忽略物质的分子结构和分子运动。我们用流动的控制方程表达基本物理守恒律的数学形式。
\begin{enumerate}
\item 流体质量守恒：
\begin{equation}\label{mass equ}
\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\vec u)=0
\end{equation}
\item 流体动量方程：
\begin{equation}\label{momen equ}
\rho \frac{\ud u}{\ud t}=\frac{\partial(-p+\tau_{xx})}{\partial x}+\frac{\partial \tau_{yx}}{\partial y}+\frac{\partial \tau_{zx}}{\partial z}
\end{equation}
\item 流体能量方程：
\begin{equation}
\label{energy equ}
\begin{split}
\rho\frac{\ud E}{\ud t}=-\nabla \cdot(p\vec u)+[\frac{\partial(u\tau_{xx})}{\partial x}+\frac{\partial(u\tau_{xy})}{\partial y}+\frac{\partial(u\tau_zz)}{\partial z}+\frac{\partial(v\tau_yx)}{\partial x}
+\frac{\partial(v\tau_yy)}{\partial y}+\\  \frac{\partial(v\tau_yz)}{\partial z}+\frac{\partial(w\tau_zx)}{\partial x}+\frac{\partial(w\tau_zy)}{\partial y}+\frac{\partial(w\tau_zz)}{\partial z}+\nabla \cdot(k\nabla T)+S_E
\end{split}
\end{equation}
\item 牛顿流体的剪应力公式:
\label{newton}
\begin{equation}
\begin{split}
\tau_{xx}=2\mu\frac{\partial u}{\partial x} +\lambda\nabla\cdot\vec u,  \tau_{xy}=\tau_{yx}=\mu(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x})\\ \dots\dots
\end{split}
\end{equation}
由以上方程\ref{mass equ},\ref{momen equ},\ref{energy equ}以及\ref{newton}可以推导出牛顿流体的Nacier-Stokes方程。对于建立有限体积法，其最有用的形式是：
\item N-s方程:
\begin{equation}
\rho \frac{\ud u}{\ud t}=-\frac{\partial p}{\partial x}+\nabla\cdot(\mu\nabla u)+S_{Mx}
\end{equation}
\begin{equation}
\rho \frac{\ud v}{\ud t}=-\frac{\partial p}{\partial y}+\nabla\cdot(\mu\nabla v)+S_{My}
\end{equation}
\begin{equation}
\rho \frac{\ud w}{\ud t}=-\frac{\partial p}{\partial z}+\nabla\cdot(\mu\nabla w)+S_{Mz}
\end{equation}
\item 对于二维不可压缩流体的流动问题，可以用流函数方程和涡量方程来描述
将流函数定义为：
\begin{displaymath}
u=\frac{\partial\psi}{\partial y},v=-\frac{\partial\psi}{\partial x}
\end{displaymath}
将涡量定义为流体质点旋转角速度的两倍，二维流动的涡量只有一个分量
\begin{displaymath}
\zeta=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}
\end{displaymath}
由二维不可压缩流体运动方程可得对流扩散方程
\begin{equation}
\frac{\partial\zeta}{\partial t}+u\frac{\partial\zeta}{\partial x}+v\frac{\partial \eta}{\partial y}=\alpha(\frac{\partial^{2}\eta}{\partial x^{2}}+\frac{\partial^{2}\eta}{\partial y^{2}})
\end{equation}
式中，$\alpha=\frac{1}{R_{e}}$为粘性扩散系数。将流函数表示的速度分量代入涡量表达式，得到流函数的泊松方程。
\begin{equation}
\frac{\partial^{2}\psi}{\partial x^{2}}+\frac{\partial^{2}\psi}{\partial y^{2}}=-\eta
\end{equation}
\end{enumerate}
\subsection{数值模拟二维不可压缩粘性流体绕圆柱的流动}
以下内容摘自\cite{cfd},是我们主要参考书目。
\includepdfmerge{cfd.pdf,198-206}
\newpage
\addcontentsline{toc}{subsection}{选录自参考书内容}
\section{程序结构说明（数据结构，主要算法，类，函数}
\subsection{后端}
\subsubsection{使用的库：UMFPACK}
UMFPACK\cite{ump}是一系列解非对称稀疏矩阵$AX=B$的库，使用了多波前算法和直接的LU分解算法。它使用ANSI/ISO C写成，提供了MATLAB前端。UMFPACK依赖于BLAS-Level 3。UMFPACK支持Windows和许多Unix平台。
解大规模线性方程组是SIMPLE算法的重点，其运算效率直接影响到整个程序的效率。后文中会提到，我们尝试了简单迭代法和松弛迭代法，但是迭代发散，所以不得不寻求别的高效解法。首先，我们的方程组中系数为零的项占绝大多数，在流函数方程组中每个方程仅有5个非零项，而在涡量方程组中不超过9个，非零项个数和方程个数n线性相关，为稀疏矩阵。对于一般矩阵，存储空间和计算时间复杂度均为$O(n^2)$,而如能利用稀疏矩阵的特性，
则可降至$O(n)$，这对于快速解出大型方程组来说是必要的。UMFPACK能充分利用矩阵的稀疏性，本身算法先进，稳定，可扩展性好，而且UMFPACK可以使用用户提供的BLAS库而非静态的算法，因此用户可以提供高效的BLAS实现（如Intel MKL）来加快运算速度。
\begin{itemize}
\item solver.h

将UMFPACK提供的C接口封装为C++类Solver，通过给定矩阵元$a_{ij}$的值和$B$便可解得X。
提供接口：
构造函数Solver(int num, int nonzero, double* right)
num为矩阵阶数，nonzero为非零项个数，right为一维数组B。
\item bool input(int i, int j, double x)

设置aij=x，如果成功返回true，如果设置非零项的总数超过了nonzero则返回false。

\item void solve()

求解方程

\item double* getSolution()

返回一维数组X。注意Solver类析构时会释放X的内存空间，所以访问时请保证Solver类未被析构

\end{itemize}
\subsubsection{主要的类}
\begin{description}
\item[\underline{cylinderDataVariant}]
cylinderDataVariant为DataVariant的一个子类，用于前端从后端获取数据。

\item[\underline{cylinderNode}]
cylinderNode用于记录每个计算节点的流函数、涡量和速度

\item[\underline{cylinderSpotStain}]
cylinderSpotStain是染色点类，用于染色点法绘制流线图。

\item[\underline{cylinderProject}]
cylinderProject是圆柱体绕流模型，提供接口如下：
\item[构造函数]cylinderProject(int l = -100, int r = 400, int u = 100, int d = -100, double dens = 0.1, double dxi = 0.1, double deta = 0.1, double dt = 0.1, double rey = 40)。
构造并初始化系统。其中l为计算坐标系的左边界，r为右边界，u为上边界，d为下边界。dens为绘制的流线密度，取值在0-1间，值越大流线越密集。dxi为计算坐标系xi方向的坐标间隔（相邻两个计算节点的横坐标差），deta为eta方向坐标间隔，dt为时间步长，rey为雷诺数。

\item[\underline{构造函数cylinderProject(const char* location)}]
读取存盘文件并继续计算。location为存盘文件路径。存盘文件可由方法dumptofile生成。

\item[\underline{DataVariant * getData(Project::DataType type, ...)}]
根据传入的数据类型传出交换数据。type可以为Project::TimeType, Project::PsiType, Project::SpotType, Project::NumberType之一，分别返回包含系统时间、流函数、染色点位置、染色点发生源数目。其中流函数需要用可变参数再提供节点的xi和eta坐标，染色点位置需要再提供染色点发生源的标号（从0到染色点发生源数目-1）。

\item[\underline{void setDensity(double dens)}void run();]
进行一步计算，系统时间增加一个时间步长

\item[\underline{void spotstainrun();}]
染色点运动一个时间步长

\item[\underline {bool dumptofile(const char* location);}]
将当前计算项目存盘，location为存盘文件路径。如果成功返回true，失败返回false。

\item[\underline{DataVariant}]
DataVariant是前端与后端进行数据交换的类型。有四种类型：Project::TimeType, Project::PsiType, Project::SpotType, Project::NumberType。提供接口如下：
\begin{itemize}
\item getX()

 仅适用于Project::PsiType, Project::SpotType
获得目标的x坐标

\item getY() \\仅适用于Project::PsiType, Project::SpotType
获得目标的y坐标

\item getZ()\\ 仅适用于Project::PsiType, Project::SpotType
获得目标的z坐标

\item getPsi()\\ 仅适用于Project::PsiType
获得节点的流函数值

\item getTime() \\仅适用于Project::TimeType
获得系统的时间

\item getNumber() \\仅适用于Project::NumberType
获得染色点发生源数目
\end{itemize}
\end{description}
\subsubsection{代码示例}
\begin{itemize}
\paragraph{简单迭代法}
\item 以下为简单迭代法算法的代码。我们之前采用这个算法，发现数值很快超过浮点数上限，算法不收敛。经过研究思考最后改用现在使用的算法。本部分代码部分采用了Intel提供的线性代数函数库\cite{intel}
\begin{lstlisting}
#pragma omp parallel for private(j, i)

    for (j = downboundary + 2; j < upboundary - 1; j++) {
        for (i = leftboundary + 2; i < rightboundary - 1; i++) {
            coordination->access(i, j).zetat = coordination->access(i, j).zeta;
        }
    }


    converge = 30;

    while (converge) {
        //TODO:converge
        #pragma omp parallel for private(j, i, node)
        for (j = downboundary + 2; j < upboundary - 1; j++) {
            for (i = leftboundary + 2; i < rightboundary - 1; i++) {
                if ((j <= 2) && (j >= -2) && (i >= leftterminal)
                    && (i <= rightterminal)) {
                    continue;
                }

                node = &coordination->access(i, j);
                node->newzetat =
                    (node->c2 * node->zeta +
                     node->c3 * (coordination->access(i, j + 2).zetat +
                                 coordination->access(i,
                                                      j + 2).zeta) +
                     node->c4 * (coordination->access(i, j + 1).zetat +
                                 coordination->access(i,
                                                      j + 1).zeta) +
                     node->c5 * (coordination->access(i, j - 1).zetat +
                                 coordination->access(i,
                                                      j - 1).zeta) +
                     node->c6 * (coordination->access(i, j - 2).zetat +
                                 coordination->access(i,
                                                      j - 2).zeta) +
                     node->c7 * (coordination->access(i + 2, j).zetat +
                                 coordination->access(i + 2,
                                                      j).zeta) +
                     node->c8 * (coordination->access(i + 1, j).zetat +
                                 coordination->access(i + 1,
                                                      j).zeta) +
                     node->c9 * (coordination->access(i - 1, j).zetat +
                                 coordination->access(i - 1,
                                                      j).zeta) +
                     node->c10 * (coordination->access(i - 2, j).zetat +
                                  coordination->access(i - 2, j).zeta)) / node->c1;
            }
        }

        j = 2;
        #pragma omp parallel for private(i, node)

        for (i = leftterminal; i <= rightterminal; i++) {
            node = &coordination->access(i, j);
            node->newzetat =
                (node->c2 * node->zeta +
                 node->c3 * (coordination->access(i, j + 2).zetat +
                             coordination->access(i,
                                                  j + 2).zeta) +
                 node->c4 * (coordination->access(i, j + 1).zetat +
                             coordination->access(i,
                                                  j + 1).zeta) +
                 node->c5 * (coordination->access(i, j - 1).zetat +
                             coordination->access(i,
                                                  j - 1).zeta) +
                 node->c6 * (cylinderBoundary->access(i, 1).zetat +
                             cylinderBoundary->access(i,
                                     1).zeta) +
                 node->c7 * (coordination->access(i + 2, j).zetat +
                             coordination->access(i + 2,
                                                  j).zeta) +
                 node->c8 * (coordination->access(i + 1, j).zetat +
                             coordination->access(i + 1,
                                                  j).zeta) +
                 node->c9 * (coordination->access(i - 1, j).zetat +
                             coordination->access(i - 1,
                                                  j).zeta) +
                 node->c10 * (coordination->access(i - 2, j).zetat +
                              coordination->access(i - 2, j).zeta)) / node->c1;
        }

        j = -2;
        #pragma omp parallel for private(i, node)

        for (i = leftterminal; i <= rightterminal; i++) {
            node = &coordination->access(i, j);
            node->newzetat =
                (node->c2 * node->zeta +
                 node->c3 * (cylinderBoundary->access(i, 0).zetat +
                             cylinderBoundary->access(i,
                                     0).zeta) +
                 node->c4 * (coordination->access(i, j + 1).zetat +
                             coordination->access(i,
                                                  j + 1).zeta) +
                 node->c5 * (coordination->access(i, j - 1).zetat +
                             coordination->access(i,
                                                  j - 1).zeta) +
                 node->c6 * (coordination->access(i, j - 2).zetat +
                             coordination->access(i,
                                                  j - 2).zeta) +
                 node->c7 * (coordination->access(i + 2, j).zetat +
                             coordination->access(i + 2,
                                                  j).zeta) +
                 node->c8 * (coordination->access(i + 1, j).zetat +
                             coordination->access(i + 1,
                                                  j).zeta) +
                 node->c9 * (coordination->access(i - 1, j).zetat +
                             coordination->access(i - 1,
                                                  j).zeta) +
                 node->c10 * (coordination->access(i - 2, j).zetat +
                              coordination->access(i - 2, j).zeta)) / node->c1;
        }

        j = 1;
        #pragma omp parallel for private(i, node)

        for (i = leftterminal; i <= rightterminal; i++) {
            node = &coordination->access(i, j);
            node->newzetat =
                    (node->c2 * node->zeta +
                     node->c3 * (coordination->access(i, j + 2).zetat +
                                 coordination->access(i,
                                                      j + 2).zeta) +
                     node->c4 * (coordination->access(i, j + 1).zetat +
                                 coordination->access(i,
                                                      j + 1).zeta) +
                     node->c5 * (cylinderBoundary->access(i, 1).zetat +
                                 cylinderBoundary->access(i, 1).zeta) +
                     node->c7 * (coordination->access(i + 2, j).zetat +
                                 coordination->access(i + 2,
                                                      j).zeta) +
                     node->c8 * (coordination->access(i + 1, j).zetat +
                                 coordination->access(i + 1,
                                                      j).zeta) +
                     node->c9 * (coordination->access(i - 1, j).zetat +
                                 coordination->access(i - 1,
                                                      j).zeta) +
                     node->c10 * (coordination->access(i - 2, j).zetat +
                                  coordination->access(i - 2, j).zeta)) / node->c1;
        }

        j = -1;
        #pragma omp parallel for private(i, node)

        for (i = leftterminal; i <= rightterminal; i++) {
            node = &coordination->access(i, j);
            node->newzetat =
                (node->c2 * node->zeta +
                 node->c4 * (cylinderBoundary->access(i, 0).zetat +
                             cylinderBoundary->access(i,
                                     0).zeta) +
                 node->c5 * (coordination->access(i, j - 1).zetat +
                             coordination->access(i,
                                                  j - 1).zeta) +
                 node->c6 * (coordination->access(i, j - 2).zetat +
                             coordination->access(i,
                                                  j - 2).zeta) +
                 node->c7 * (coordination->access(i + 2, j).zetat +
                             coordination->access(i + 2,
                                                  j).zeta) +
                 node->c8 * (coordination->access(i + 1, j).zetat +
                             coordination->access(i + 1,
                                                  j).zeta) +
                 node->c9 * (coordination->access(i - 1, j).zetat +
                             coordination->access(i - 1,
                                                  j).zeta) +
                 node->c10 * (coordination->access(i - 2, j).zetat +
                              coordination->access(i - 2, j).zeta)) / node->c1;
        }


        #pragma omp parallel for private(j, i)

        for (j = downboundary + 2; j < upboundary - 1; j++) {
            for (i = leftboundary + 2; i < rightboundary - 1; i++) {
                coordination->access(i, j).zetat =
                    coordination->access(i, j).newzetat;
            }
        }

        converge -= 1;
    }
    */

    /* flush zetat back to zeta 
    #pragma omp parallel for private(j, i)

    for (j = downboundary + 2; j < upboundary - 1; j++) {
        for (i = leftboundary + 2; i < rightboundary - 1; i++) {
            coordination->access(i, j).zeta = coordination->access(i, j).zetat;
        }
    } */
\end{lstlisting}
\paragraph{多波前算法}
\item 以下列出多波前算法计算$\zeta$的代码，是我们最终采用的版本。它用来求解大规模稀疏矩阵，在我们的课题里用来求解各个点的$eta$值。在经历前一种算法的挫折后，我们很惊喜地发现在当前算法下数据很快就收敛了。
\begin{lstlisting}
void cylinderProject::calculateNewZeta()
{
    double vxi, veta;
    double uxi, ueta;
    double hxi, heta;
    double lambdaxi, lambdaeta;
    int i, j;
    cylinderNode * node;
    /* Calculating new zeta at t + deltat on inner nodes */
    /* Calculating coefficients */
    #pragma omp parallel for private(i, j, node, hxi, heta, vxi, veta, uxi, ueta, lambdaxi, lambdaeta)

    for (i = downboundary + 2; i < upboundary - 1; i++) {
        for (j = leftboundary + 2; j < rightboundary - 1; j++) {
            if ((i <= 1) && (i >= -1) && (j >= leftterminal)
                && (j <= rightterminal)) {
                /* Near or on the cylinder, do nothing here */
                continue;
            }

            node = &coordination->access(j, i);
            hxi = node->hxi;
            heta = node->heta;
            vxi =
                (coordination->access(j, i + 1).psi -
                 coordination->access(j, i - 1).psi) / (2 * heta * deltaeta);
            veta =
                -(coordination->access(j + 1, i).psi -
                  coordination->access(j - 1, i).psi) / (2 * hxi * deltaxi);
            uxi = vxi / hxi;
            ueta = veta / heta;

            if (Re < 1000) {
                lambdaxi = 1 - 1 / exp(abs(uxi));
                lambdaeta = 1 - 1 / exp(abs(ueta));
            } else {
                lambdaxi = 1;
                lambdaeta = 1;
            }

            node->c1 =
                -2 / deltat - lambdaxi * abs(uxi) / (2 * deltaxi) -
                lambdaeta * abs(ueta) / (2 * deltaeta) -
                4 / (hxi * hxi * deltaxi * deltaxi * Re) -
                4 / (heta * heta * deltaeta * deltaeta * Re);
            node->c2 =
                -2 / deltat + lambdaxi * abs(uxi) / (2 * deltaxi) +
                lambdaeta * abs(ueta) / (2 * deltaeta) +
                4 / (hxi * hxi * deltaxi * deltaxi * Re) +
                4 / (heta * heta * deltaeta * deltaeta * Re);
            node->c3 = -(ueta - lambdaeta * abs(ueta)) / (12 * deltaeta);
            node->c4 =
                (2 * ueta - lambdaeta * abs(ueta)) / (3 * deltaeta) -
                2 / (heta * heta * deltaeta * deltaeta * Re);
            node->c5 =
                -(2 * ueta + lambdaeta * abs(ueta)) / (3 * deltaeta) -
                2 / (heta * heta * deltaeta * deltaeta * Re);
            node->c6 = (ueta + lambdaeta * abs(ueta)) / (12 * deltaeta);
            node->c7 = -(uxi - lambdaxi * abs(uxi)) / (12 * deltaxi);
            node->c8 =
                (2 * uxi - lambdaxi * abs(uxi)) / (3 * deltaxi) -
                2 / (hxi * hxi * deltaxi * deltaxi * Re);
            node->c9 =
                -(2 * uxi + lambdaxi * abs(uxi)) / (3 * deltaxi) -
                2 / (hxi * hxi * deltaxi * deltaxi * Re);
            node->c10 = (uxi + lambdaxi * abs(uxi)) / (12 * deltaxi);
        }
    }

    #pragma omp parallel for private(j, hxi, heta, vxi, veta, uxi, ueta, lambdaxi, lambdaeta, node)

    for (j = leftterminal; j <= rightterminal; j++) {
        /* Upper half */
        node = &coordination->access(j, 1);
        hxi = node->hxi;
        heta = node->heta;
        vxi =
            (coordination->access(j, 2).psi -
             cylinderBoundary->access(j, 1).psi) / (2 * heta * deltaeta);
        veta =
            -(coordination->access(j + 1, 1).psi -
              coordination->access(j - 1, 1).psi) / (2 * hxi * deltaxi);
        uxi = vxi / hxi;
        ueta = veta / heta;

        if (Re < 1000) {
            lambdaxi = 1 - 1 / exp(abs(uxi));
            lambdaeta = 1 - 1 / exp(abs(ueta));
        } else {
            lambdaxi = 1;
            lambdaeta = 1;
        }

        node->c1 =
            -2 / deltat - lambdaxi * abs(uxi) / (2 * deltaxi) -
            (-ueta) / (2 * deltaeta) -
            4 / (hxi * hxi * deltaxi * deltaxi * Re) -
            4 / (heta * heta * deltaeta * deltaeta * Re);
        node->c2 =
            -2 / deltat + lambdaxi * abs(uxi) / (2 * deltaxi) +
            (-ueta) / (2 * deltaeta) +
            4 / (hxi * hxi * deltaxi * deltaxi * Re) +
            4 / (heta * heta * deltaeta * deltaeta * Re);
        node->c3 = -(ueta - (-ueta)) / (12 * deltaeta);
        node->c4 =
            (2 * ueta - (-ueta)) / (3 * deltaeta) -
            2 / (heta * heta * deltaeta * deltaeta * Re);
        node->c5 =
            -(2 * ueta + (-ueta)) / (3 * deltaeta) -
            2 / (heta * heta * deltaeta * deltaeta * Re);
        node->c6 = (ueta + (-ueta)) / (12 * deltaeta);
        node->c7 = -(uxi - lambdaxi * abs(uxi)) / (12 * deltaxi);
        node->c8 =
            (2 * uxi - lambdaxi * abs(uxi)) / (3 * deltaxi) -
            2 / (hxi * hxi * deltaxi * deltaxi * Re);
        node->c9 =
            -(2 * uxi + lambdaxi * abs(uxi)) / (3 * deltaxi) -
            2 / (hxi * hxi * deltaxi * deltaxi * Re);
        node->c10 = (uxi + lambdaxi * abs(uxi)) / (12 * deltaxi);

        /* Lower half */
        node = &coordination->access(j, -1);
        hxi = node->hxi;
        heta = node->heta;
        vxi =
            (cylinderBoundary->access(j, 0).psi -
             coordination->access(j, -2).psi) / (2 * heta * deltaeta);
        veta =
            -(coordination->access(j + 1, -1).psi -
              coordination->access(j - 1, -1).psi) / (2 * hxi * deltaxi);
        uxi = vxi / hxi;
        ueta = veta / heta;

        if (Re < 1000) {
            lambdaxi = 1 - 1 / exp(abs(uxi));
            lambdaeta = 1 - 1 / exp(abs(ueta));
        } else {
            lambdaxi = 1;
            lambdaeta = 1;
        }

        node->c1 =
            -2 / deltat - lambdaxi * abs(uxi) / (2 * deltaxi) -
            ueta / (2 * deltaeta) -
            4 / (hxi * hxi * deltaxi * deltaxi * Re) -
            4 / (heta * heta * deltaeta * deltaeta * Re);
        node->c2 =
            -2 / deltat + lambdaxi * abs(uxi) / (2 * deltaxi) +
            ueta / (2 * deltaeta) +
            4 / (hxi * hxi * deltaxi * deltaxi * Re) +
            4 / (heta * heta * deltaeta * deltaeta * Re);
        node->c3 = -(ueta - ueta) / (12 * deltaeta);
        node->c4 =
            (2 * ueta - ueta) / (3 * deltaeta) -
            2 / (heta * heta * deltaeta * deltaeta * Re);
        node->c5 =
            -(2 * ueta + ueta) / (3 * deltaeta) -
            2 / (heta * heta * deltaeta * deltaeta * Re);
        node->c6 = (ueta + ueta) / (12 * deltaeta);
        node->c7 = -(uxi - lambdaxi * abs(uxi)) / (12 * deltaxi);
        node->c8 =
            (2 * uxi - lambdaxi * abs(uxi)) / (3 * deltaxi) -
            2 / (hxi * hxi * deltaxi * deltaxi * Re);
        node->c9 =
            -(2 * uxi + lambdaxi * abs(uxi)) / (3 * deltaxi) -
            2 / (hxi * hxi * deltaxi * deltaxi * Re);
        node->c10 = (uxi + lambdaxi * abs(uxi)) / (12 * deltaxi);
    }
\end{lstlisting}
\end{itemize}


\subsection{前端}

\subsubsection{使用的库 ：Qt}
Qt\cite{qt}是一个跨平台的C++应用程序开发框架。广泛用于开发GUI程序，这种情况下又
被称为部件工具箱。也可用于开发非GUI程序，比如控制台工具和服务器。

Qt是自由且开放源代码的软件，在GNU较宽松公共许可证条款下发布。所有版本
都支持广泛的编译器，包括GCC的C++编译器和Visual Studio。

\subsubsection{主要的类}

在前端程序中较为重要的类有：

\begin{description}
\item[\underline{ProjectMainWindow}]
  ProjectMainWindow是Qt类QMainWindow的子类，是图形界面前端的主要组成部
  分，即主窗口对象。
\item[\underline{CalThread}]
  CalThread类是Qt类QThread的子类，是对后端计算进行Qt包装的主类。为了让
  后端运算不影响前端的访问响应，在这个类里把运算封装到Qt的线程类
  QThread中，让计算过程在主线程以外的线程运行。
\item[\underline{DisplayWidget}]
  DisplayWidget类是Qt类QGraphicsView的子类，是负责将后端生成的数据图像
  化的主类。这个类使用了OpenGL引擎来渲染作图，也同时负责将后端生成的数
  据缓存起来用于前端的反复访问。
\item[\underline{InputWidget}]
  InputWidget类是一个GUI组件，主要用于收集计算所需要的初始条件并负责把
  这些初始条件传递给后端用以初始化。
\item[\underline{ControlWidget}]
  ControlWidget类是一个GUI组件，主要用于控制图像显示中的相关参数和控制
  后端计算的具体流程。
\item[\underline{SudokuWidget}]
  SudokuWidget是Qt类QWidget的子类，是数独游戏的主窗口。
\end{description}

\subsection{项目整体组织}

\subsubsection{项目构建管理}

本项目使用qmake作为项目的构建管理工具。

qmake\cite{qt}是一个协助简化跨平台进行专案开发的构建过程的工具程序，Qt附带的工
具之一 。qmake能够自动生成Makefile、Microsoft Visual Studio 项目文件
和 xcode 项目文件。不管源代码是否是用Qt写的，都能使用qmake，因此qmake
能用于很多软件的构建过程。

本项目使用qmake的考虑主要是由于其对Qt的支持是内建的，因此在兼容性方面它应
当有最好的表现。

\subsubsection{项目结构管理}
本项目采用了模块化的组织方式。

项目的主体代码在文件夹src/下，组织在各个不同的文件夹里，代表不同的模块。
每个模块在编译时各自生成静态库文件，保存在lib/下，最后编译main模块时，
将各模块部分链接在一起生成最终的可执行文件并保存在bin/下。

这样的组织结构使得程序高度灵活化，每一个模块都可以单独拿出作为一个库使
用，而所有的模块作为库的对外借口文件都保存在include/下。

doc/下保存了项目的全部文档文件，demo/下是项目早期用于验证算法思路的一
个实验性程序。

src/下的模块有：
\begin{itemize}
\item 2D模块，用于计算圆柱体绕流的数据，是后端的主要组成模块。
\item frontend模块，用于图形化2D模块计算出来的数据，是前端的主要组成模
  块。
\item sudoku模块，这是数独游戏的主要后端构成模块，是用纯C语言写的。
\item sudokugui模块，这是数独游戏的Qt GUI前端模块。
\end{itemize}

\subsubsection{项目数据管理}

本项目使用git\cite{git}来完成项目的数据管理与版本管理工作，并使用
googlecode的服务共享代码实现分布式开发。

Git是一个由Linux创始人，Linus Torvalds为了更好地管理linux内核开发而创
立的分布式版本控制／软件配置管理软件。


\section{课题过程及结果分析}
\subsection{课题完成情况}
在开学第二次计算概论课上，我们几个人就觉得组成小组，完成课题。由于基础较好，对课题的前景十分乐观。我们挑选了很久项目的主题，最终决定做计算流体力学相关的课题。一来比较有挑战性，二来也切合物理学科。我们信心百倍的完成了开题报告，但在接下来的时间里只是看参考资料，相关书目。真正开始做已经是11月初。由于课题难度之大，直到近日才初步做完，可谓蹉跎。

具体分工来看，徐智怡完成后端的大部分工作。沈钟灵负责数学物理的相关问题、部分后端程序以及项目报告的书写。黄康靖同学则完成了前端的工作。在这个项目上，三人都投入大量时间精力，同时也收获了很多。

最后的效果差强人意，远不及当初那么理想。甚至还存在一些我们暂时无法解决
的bug.比如圆柱的大小显得很不合适。以下将列出在完成项目过程中我们所遇到的问题，以及项目的一些优点特点。
\subsection{项目的问题和优点}
\subsubsection{问题}
\begin{itemize}
\item 使用迭代法计算流函数和涡量时不收敛，应该是算法的精度不够高而至。最后我们换了方程的迭代方法，并利用umfpack求解。得到了较好的结果。
\item 不能从数值上判断收敛，只能从图像上看。在作出前端以后我们才对初始值设定进行了进一步的修正，使得结果有更高可信度。
\item 由于直接调用OpenGL接口做图过于繁琐冗长且难以调试，碍于时间有限的
  缘故，我们使用了Qt的直接接口辅以Qt自带的OpenGL引擎渲染完成了作图，
  并且放弃了原先作等流函数图的计划，只作出了流线图。

\item 使用迭代法计算流函数和涡量时不收敛，最后用 多波前迭代，使用umfpack库解决

\item 不能从数值上判断收敛，只能从图像上看。所以最终对数据又有一些修正。


\end{itemize}
\subsubsection{亮点}
\begin{itemize}
\item 由于参考资料内容不够详尽，很多算法、数理方面的推导由自己完成。例如：圆柱体边界特殊处理方法，上下边界左右节点比较特殊，不是正交的。我们认为它是连续，而涡量则是边界的上下极限平均值。尾流约束条件没有用微分约束，认为是均匀量。又如，所采用的保角变换的拉梅系数的推导等等。
\item 使用了Git作为项目的版本控制工具，便于检索和比较历史版本，有利于
  增加开发效率。
\item 使用了Openmp进行并行计算。OpenMp是由OpenMP Architecture Review Board牵头提出的，并已被广泛接受的，用于共享内存并行系统的多线程程序设计的一套指导性注释(Compiler Directive)。OpenMP支持的编程语言包括C语言、C++和Fortran；而支持OpenMp的编译器包括Sun Compiler，GNU Compiler和Intel Compiler等。OpenMp提供了对并行算法的高层的抽象描述，程序员通过在源代码中加入专用的pragma来指明自己的意图，由此编译器可以自动将程序进行并行化，并在必要之处加入同步互斥以及通信。当选择忽略这些pragma，或者编译器不支持OpenMp时，程序又可退化为通常的程序(一般为串行)，代码仍然可以正常运作，只是不能利用多线程来加速程序执行。
\item 可以中途存盘(仅后端可以，前端暂未实现）

\end{itemize}


\begin{thebibliography}
{99}
 \bibitem{cfd} 李万平
 \emph{计算流体力学},华中科技大学 
(2004)
\bibitem{intel}Intel Math Kernal Library
\bibitem{ump}UMPACK User Guide
\bibitem{qt}Qt 4.8 reference Documentation
\bibitem{git}Git User Manual


  \end{thebibliography}




\end{document}  